Double scatter simulation for improved reconstruction of positron emission tomography data

ABSTRACT

Methods for simulating, and correcting for, doubly scattered annihilation gamma-ray photons in both time-of-flight (TOF) and non-TOF positron emission tomography scan data are disclosed.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority under 35 U.S.C. § 119(e) to U.S.Provisional Application No. 62/757,893, filed Nov. 9, 2018, the entirecontents of which are incorporated herein by reference.

FIELD

The present disclosure generally relates to nuclear medicine, andsystems for obtaining nuclear medicine images. In particular, thepresent disclosure relates to systems and methods for reconstructingnuclear medicine images from time-of-flight (TOF) positron emissiontomography (PET) data.

BACKGROUND

Nuclear medicine is a unique medical specialty wherein radiation is usedto acquire images which show the function and anatomy of organs, bonesor tissues of the body. Radiopharmaceuticals are introduced into thebody, either by injection or ingestion, and are attracted to specificorgans, bones or tissues of interest. Such radiopharmaceuticals producegamma photon emissions which emanate from the body and are captured by ascintillation crystal, with which the photons interact to produceflashes of light or “events.” Events are detected by an array ofphotodetectors, such as photomultiplier tubes, and their spatiallocations or positions are calculated and stored. In this way, an imageof the organ or tissue under study is created from detection of thedistribution of the radioisotopes in the body.

One particular nuclear medicine imaging technique is known as PositronEmission Tomography, or PET. PET is used to produce images fordiagnosing the biochemistry or physiology of a specific organ, tumor orother metabolically active site. Measurement of the tissue concentrationof a positron emitting radionuclide is based on coincidence detection ofthe two gamma photons arising from positron annihilation. When apositron is annihilated by an electron, two 511 keV gamma photons aresimultaneously produced and travel in approximately opposite directions.Gamma photons produced by an annihilation event can be detected by apair of oppositely disposed radiation detectors capable of producing asignal in response to the interaction of the gamma photons with ascintillation crystal. Annihilation events are typically identified by atime coincidence between the detection of the two 511 keV gamma photonsin the two oppositely disposed detectors, i.e., the gamma photonemissions are detected virtually simultaneously by each detector. Whentwo oppositely disposed gamma photons each strike an oppositely disposeddetector to produce a time coincidence event, they also identify a lineof response, or LOR, along which the annihilation event has occurred. Anexample of a PET method and apparatus is described in U.S. Pat. No.6,858,847, which patent is incorporated herein by reference in itsentirety.

After being integrated and sorted into LORs defined by the positions ofthe detectors in the PET camera, the coincidence event data are used toreconstruct a three-dimensional distribution of the positron-emittingradionuclide within the patient. In two-dimensional PET, each 2Dtransverse (perpendicular to the axis of the PET scanner) section or“slice” of the radionuclide distribution is reconstructed independentlyof adjacent sections, using only LORs that are (approximately)perpendicular to the axis of the scanner (referred to as the z axis). Infully three-dimensional PET, nearly all the LOR data are used for thereconstruction. The positions of these LORs in space are characterizedby their radial distance, s, from the z axis, their azimuthal angle, Φ,around the z axis, their polar angle, Θ, with respect to the z axis, andthe z position of their closest approach to the z axis. These LOR dataare typically arranged into a set of “sinograms”, p(s, Φ; Θ, z), which,for fixed values of Θ and z, represents a two dimensional parallelprojection of the three dimensional radionuclide distribution within thepatient. All of the LORs in the sinogram p(s, Φ) having fixed values ofΘ and z are essentially co-planar. In this format, a single fixed pointin the emitter distribution f(x,y,z) that falls within this (Θ, z) planetraces a sinusoid in the sinogram. In each sinogram, there is one rowcontaining the LORs for a particular azimuthal angle Φ; each such rowcorresponds to a one-dimensional parallel projection of the tracerdistribution at a different projection angle. This is shown conceptuallyin FIG. 1 for case in which the plane of the sinogram is perpendicularto the z axis.

An event is registered if both crystals detect an annihilation photonwithin a coincidence time window τ (e.g., on the order of 4-5 ns),depending on the timing properties of the scintillator and the field ofview. Aside from the effect of photon scatter, as discussed below, apair of detectors is sensitive only to coincidence events originating inthe volume between the two detectors, thereby eliminating the need forphysical collimation, and thus significantly increasing sensitivity.Accurate corrections can be made for the self-absorption of photonswithin the patient (i.e., attenuation correction) so that accuratemeasurements of tracer concentration can be made.

The number of time coincidences detected per second within a field ofview (FOV) of a detector is the count rate of the detector. The countrate at each of two oppositely disposed detectors, A and B, can bereferred to as singles counts, or singles, S_(A) and S_(B). The timerequired for a gamma photon to travel from its point of origin to apoint of detection is referred to as the time of flight, or TOF, of thegamma photon. TOF is dependent upon the speed of light c and thedistance traveled. A time coincidence, or coincidence event, isidentified if the time difference between the arrival of signals in apair of oppositely disposed detectors is within the coincidence timewindow τ. In conventional PET, the coincidence detection time window τis wide enough so that an annihilation event occurring anywhere withinthe object would produce annihilation gamma photons reaching theirrespective detectors within the coincidence window. Coincidence timewindows of 4.5-12 nsec are common for conventional PET, and are largelydetermined by the time resolution capabilities of the detectors andelectronics.

As illustrated in FIG. 2, if an annihilation event occurs at themidpoint of a LOR, the TOF of the gamma photon detected in detector A(T_(A)) is equal to the TOF of the gamma photon detected in detector B(T_(B)). If an annihilation event occurs at a distance Δx from themidpoint of the LOR, the difference between T_(A) and T_(B) is Δt=2Δx/c,where c is the speed of light. If d is the distance between thedetectors, the TOF difference Δt could take any value from −d/c to +d/c,depending on the location of the annihilation event.

In contrast to conventional PET, TOF-PET is based on recording thedifference Δt between the detection times of the two gamma photonsarising from the positron annihilation event in sub-intervals of thetotal coincidence window τ. This measurement allows the annihilationevent to be localized along the LOR with a resolution of about 30-180 mmFWHM, assuming a time resolution of 200-1200 ps (picoseconds). Thoughless accurate than the spatial resolution of the scanner, thisapproximate localization is effective in reducing noise contributionsboth from random coincidence events and from scattered and unscatteredphoton coincidences that actually originated elsewhere in the object.This improves both the stability of the reconstruction and thesignal-to-noise ratio (SNR) in the final image, especially when imaginglarge objects. TOF-PET may increase image SNR by a factor of 2 or morecompared to conventional PET.

TOF scanners developed in the early 1980s were used for research andclinical applications, but the SNR gain provided by the TOF measurementsof about 500 ps resolution was offset by poorer spatial resolution andlower sensitivity due to the low stopping power of the BaF₂ and CsFscintillation crystals used in such scanners. Consequently, those TOFsystems could not compete successfully with conventional (non-TOF) BGOscanners. As a result, TOF-PET almost completely disappeared from thescene in the 1990s. Today, faster electronics and crystals such as LSOand LaBr₃ reopen the prospect of exploiting the TOF information withoutcompromising other parameters such as the count rate, the sensitivity,and the energy and spatial resolutions.

Septumless or “3D” PET scanners (i.e., without interplane septa)currently constitute a large percentage of the total market for PETimaging. Because of the lack of interplane septa, scattered events(i.e., annihilation photons undergoing Compton scattering beforereaching the detector) may represent a large portion of the measureddata (e.g., up to 50% or more in clinical studies). An example of suchscatter is shown in FIG. 3, which illustrates a septumless PET scannerutilizing a ring detector configuration, which is also applicable foruse with the present invention.

An annihilation event occurring at emission point 21 produces twooppositely traveling gamma photons along LOR 22. One of the gammaphotons, however, may undergo Compton scattering at scatter point 23,which changes its travel direction to path 24. Consequently, while thefirst gamma photon is detected by detector A in line with theoriginating LOR 22, the scattered gamma photon will be detected indetector B, rather than by detector C. Consequently, the coincidenceevent detected in detectors A and C will result in a false LOR 25 beingidentified, instead of the correct LOR 22.

An example of techniques to correct for such scattering in TOF PETutilizing a single scatter simulation (SSS) algorithm is disclosed inU.S. Pat. No. 7,397,035, the contents of which are incorporated hereinby reference in entirety.

However, in reality, the scattering of gamma photons from annihilationevents before reaching the detectors also includes multiple Comptonscattering. For modern scanners with good energy resolution and a narrowphotopeak energy window, recent Monte Carlo simulations have shown thatonly double scatter makes a significant contribution to the multiplescatter and higher order scatter contributions are very small ornegligible.

Currently, only the single scatter component is accurately modeled, andmultiple scatter is accounted for by scaling this model of singlescatter to the measured data. This is problematic for several reasons,including the possibility that single and multiple scatter aredistributed differently. Additionally, the regions of the dataappropriate for use in such scaling can be difficult to identify, andthe measured data in those regions are often too sparse or noisy tosupport an accurate scaling.

Therefore, there is a need for an improved scatter simulation method toaccurately model single plus double scatter contributions to TOF-PETdata for more accurate approximation of the total scatter contributions.

SUMMARY

The methods disclosed herein provide a physically complete, accurate,and absolutely scaled model for double scatter estimation in a PETscanner. Also disclosed herein is a novel sampling scheme for the doublescatter integration that improves the efficiency of the calculation.

According to an aspect of the present disclosure, a computer-implementedmethod for applying scatter correction to a scan data acquired in a TOFPET scanner is disclosed. The method comprises: (A) selecting a pair ofdetector positions in the TOF PET scanner's detector ring; (B)numerically estimating the double scatter coincidence event rate ofannihilation photons in the selected pair of detector positions, whereinthe annihilation photons in a double scatter coincidence event areassumed to be scattered once, each photon by a different scatter pointin the TOF PET scanner's image volume, during a TOF PET acquisition scanperiod; (C) selecting another pair of detector positions; (D) repeatingstep (B) wherein the selected pair of detector positions is the anotherpair of detector positions; (E) repeating steps (C) and (D) until alldetector pairs of interest have been selected; and (F) scattercorrecting acquired scan data obtained from the pair of detectorpositions during the TOF PET acquisition scan period to reconstruct aTOF PET image based on the scan data.

According to some embodiments, a method for applying scatter correctionto a scan data acquired in a TOF PET scanner comprises: (A) selecting apair of detector positions in the TOF PET scanner's detector ring; (B)numerically estimating the double scatter coincidence event rate ofannihilation photons in the selected pair of detector positions, whereinone of the pair of annihilation photons in a double scatter coincidenceevent is assumed not to be scattered and the other of the two photons isscattered twice by a scatter point in the TOF PET scanner's imagevolume, during a TOF PET acquisition scan period; (C) selecting anotherpair of detector positions; (D) repeating step (B) wherein the pair ofdetector positions is the another pair of detector positions; and (E)repeating steps (C) and (D) until all detector pairs of interest havebeen selected; and (F) scatter correcting the acquired scan dataobtained from the pair of detector positions during the TOF PETacquisition scan period to reconstruct a TOF PET image based on the scandata.

According to some embodiments, a system for processing andreconstructing TOF PET sinogram data is disclosed. The system comprises:a processor capable of executing instructions; and a non-transitory,machine readable storage medium encoded with program instructions forcontrolling a TOF PET scanner, such that when the processor executes theprogram instructions, the processor performs at least one of thedisclosed methods for applying scatter correction to a scan dataacquired in the TOF PET scanner.

Also disclosed is a non-transitory, machine readable storage mediumencoded with program instructions for controlling a TOF PET scanner,such that when a processor executes the program instructions, theprocessor performs at least one of the disclosed methods for applyingscatter correction to a scan data acquired in the TOF PET scanner.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will now be more fully described by way of example withreference to the accompanying drawings in which:

FIG. 1 is a diagram illustrating the relationship between PET projectiondata and a sinogram;

FIG. 2 is a diagram illustrating the concept of time of flight in PETimaging;

FIG. 3 is a diagram illustrating an example of a trajectory for a singlescattering coincidence event.

FIG. 4 is a diagram illustrating a possible trajectory for doublescattering coincidence event with the three possible double scattercontributions to LOR (D₁, D₂) depending on the location of the emissionsource on the trajectory.

FIG. 5 is a diagram illustrating a possible first scatter point S₁ andthree possible second scatter points, S₂, S′₂ and S″₂, defining 3possible double scatter pair samples.

FIG. 6 is a diagram illustrating a supercell grid for N=16 and M=4.

FIG. 7 is a diagram illustrating the 8 double scatter pairs associatedwith one supercell.

FIG. 8 is a plot of average sinogram radial profiles comparingvalidation data against Monte Carlo simulations for a 20 cm cylinderphantom.

FIG. 9 is a plot of average sinogram radial profiles comparingvalidation data against Monte Carlo simulations for a 30 cm cylinderphantom.

FIGS. 10A-10D are coronal images from four different reconstructions ofclinical PET data using the radiotracer ¹⁸F-DCFPyL. ‘SSS’=single scattersimulation, ‘DSS’=double scatter simulation, ‘Rel’=relative scaling,‘Abs’=absolute scaling. ‘Relative total activity’ refers to the totalactivity in the imaged portion of the body compared to the ‘SSS+DSS Abs’case.

FIGS. 11A-11C are coronal images from three different reconstructions ofclinical PET data using the radiotracer ¹⁸F-DCFPyL (same patient as inFIGS. 10A-10D). ‘SF’=scatter fraction, ‘it’=number of iterations used inthe reconstruction, ‘sec’=seconds required for the reconstruction,‘dbl_mult’=multiplicity factor used for selecting double scatter points.

FIGS. 12A-12D are coronal images from four different reconstructions ofa second clinical PET study using the radiotracer ¹⁸F-DCFPyL. Thenotations have the same meaning as in FIGS. 10A-10D.

FIGS. 13A-13C are coronal images from three different reconstructions ofa clinical PET study using the radiotracer ¹⁸F-FDG. The notations havethe same meaning as in FIGS. 10A-10D. ‘M’=the multiplicity factor usedfor selecting double scatter points.

FIGS. 14A-14D are coronal images from two different reconstructions foreach of two clinical PET studies using the radiotracer ¹⁸F-FDG. Thenotations have the same meaning as in FIGS. 10A-10D.

FIGS. 15A-15D are coronal images from two different reconstructions foreach of two clinical PET studies using the radiotracer ⁶⁸Ga-PSMA. Thenotations have the same meaning as in FIGS. 10A-10D.

FIGS. 16A-16B are coronal images from two different reconstructions of aclinical PET study using the radiotracer ⁶⁸Ga-PSMA. The notations havethe same meaning as in FIGS. 10A-10D.

FIG. 17 is a schematic illustration of a transaxial slice of thesimulated phantom.

FIG. 18A shows a comparison of Monte Carlo simulation double scatter TOFsinograms to the corresponding DSS TOF sinograms for 8 TOF bins withsignificant data.

FIG. 18B show a comparison of Monte Carlo simulation total scatter (allorders) TOF sinograms to the corresponding SSS+DSS TOF sinograms.

FIG. 19 is a high-level block diagram of a PET scanner system.

DETAILED DESCRIPTION

This description of the exemplary embodiments is intended to be read inconnection with the accompanying drawings, which are to be consideredpart of the entire written description. The disclosed embodiments aremerely exemplary of the invention and the invention may be embodied invarious and alternative forms.

Disclosed herein is a new double scatter simulation (DSS) algorithmthat, when combined with the previously known single scatter simulation(SSS) technique, accurately approximates the total Compton scatteringcontribution to PET data for both TOF and non-TOF PET scanning. Thesedouble scatter contributions are efficiently computed by considering asubset of pairs of the single scatter points. By fully accounting forthe physics, an absolute scaling is achieved, reducing or eliminatingthe need for a sometimes problematic scaling to the emission data.

The DSS algorithm disclosed herein is predicated on the prior existenceof an initial PET emission image (positron emissions/volume/time) of theperson or object being scanned in the PET camera, and the correspondingimage of the linear attenuation coefficients for the annihilationradiation within that person or object. This initial emission image maybe uncorrected for scattered radiation, and thus inaccurate. Thecontribution of scattered radiation to this image, estimated by means ofthe SSS+DSS simulations, can be used to more accurately model theemission data, and thereby allow a more accurate reconstruction of theemission image. This process of scatter estimation followed by acorrected reconstruction can be iterated until the scattered radiationis fully accounted for.

FIG. 4 shows a schematic illustration of a possible trajectory for thegamma photons from a double scattering coincidence event. D₁ and D₂ aredetectors on a detector ring. S₁ and S₂ are possible scatter points.There are three distinct contributions, depending on the location of theemission point (the gamma annihilation point) relative to the twoscatter points. For emission points on the rays (S₁, D₁) and (S₂, D₂)one photon will be unscattered and the other will be scattered twice.For emission points on the ray (S₁, S₂) both photons will be scatteredonce. According to an embodiment of the present disclosure, the totalcontribution from double scattering to the LOR (D₁, D₂) in a TOF offsetbin n, from emission points along the ray R¹¹ between S₁ and D₁, isestimated as in the following equation (1):

$\begin{matrix}\begin{matrix}{R_{n}^{11} = {\int\limits_{V_{S}}{\int\limits_{V_{S}}{\left\lbrack {\int_{D_{1}}^{S_{1}}{{h_{n}\left\lbrack {\Delta \; {t(s)}} \right\rbrack}{n_{e}(s)}\; {ds}}} \right\rbrack_{1}\left\lbrack \frac{\sigma_{D_{1}}{dV}_{S_{1}}^{2/3}}{2\pi \; R_{S_{1}D_{1}}^{2}} \right\rbrack}_{2}}}} \\{{\left\lbrack e^{\int_{D_{1}}^{S_{1}}{{\mu {(E_{0,s})}}{ds}}} \right\rbrack_{3}\left\lbrack {{\mu \left( {E_{0},S_{1}} \right)}{dV}_{S_{1}}^{1/3}} \right\rbrack}_{4}} \\{{\left\lbrack {\frac{1}{\sigma_{c}\left( E_{0} \right)}\frac{d\; {\sigma_{c}\left( {E_{0},E_{1}} \right)}}{d\; \Omega}} \right\rbrack_{5}\left\lbrack \frac{{dV}_{S_{2}}^{2/3}}{R_{S_{1}S_{2}}^{2}} \right\rbrack}_{6}} \\{{\left\lbrack e^{- {\int_{S_{1}}^{S_{2}}{{\mu {({E_{1},s})}}{ds}}}} \right\rbrack_{7}\left\lbrack {{\mu \left( {E_{1},S_{2}} \right)}{dV}_{S_{2}}^{1/3}} \right\rbrack}_{8}} \\{{\left\lbrack {\frac{1}{\sigma_{c}\left( E_{1} \right)}\frac{d\; {\sigma_{c}\left( {E_{1},E_{12}} \right)}}{d\; \Omega}} \right\rbrack_{9}\left\lbrack \frac{\sigma_{D_{2}}}{R_{S_{2}D_{2}}^{2}} \right\rbrack}_{10}} \\{{\left\lbrack e^{- {\int_{S_{2}}^{D_{2}}{{\mu {({E_{12},s})}}{ds}}}} \right\rbrack_{11}\left\lbrack {{\epsilon_{D_{1}}\left( E_{0} \right)}{\epsilon_{D_{2}}\left( E_{12} \right)}} \right\rbrack}_{12}}\end{matrix} & (1)\end{matrix}$

Only Compton scattering from free electrons, with total cross-sectionσ_(c) and differential scattering cross-section dσ_(c)/dΩ, is modeled.E₀ is the initial annihilation photon energy (511 keV). If θ_(i) is thescattering angle at S_(i), the energy of a singly scattered photon isE_(i)=E₀/(2−cos θ_(i)) and the energy of a doubly scattered photon isE₁₂=E₀/(3−cos θ₁−cos θ₂). E₁₂ does not depend on the order of the twoscatter events. There is a double volume integral over the two scatterpoints, and dV_(S) is the incremental scattering volume. n_(e)(s) is theemission rate density (annihilations/volume/time) at the spatialposition s and μ(E, s) is the linear attenuation coefficient for energyE at s. R_(SD) and R_(SS) are the detector-scatter and scatter-scatterray lengths, respectively. ϵ_(D)(E) is the detection efficiency atenergy E, and σ_(D) is the geometrical cross section of the detector.Δt(s) is the coincidence detection time offset for two photons emittedfrom the points, which follow the double scatter trajectory. It is equalto their path length difference divided by the speed of light. h_(n)[ ]is an index function whose value is 1 if Δt(s) falls in the n^(th)discrete TOF offset bin, and 0 otherwise. The time resolution of thedetectors is not modeled at this stage, but is applied after allcontributions to the LOR (D₁, D₂) have been computed. The number andwidths of the discrete TOF offset bins may be chosen so that they fullycover the total coincidence window τ. For the non-TOF case, there isonly a single time bin to consider, which is the total coincidencewindow τ, and the function h_(n)[ ] is always 1, so the time offsetsΔt(s) do not need to be computed.

Referring to the numbered brackets in equation (1), each numberedbracket being referred to as a term in the following discussion, Term 1is the integral of emission density along (D₁, S₁). Term 2 is thegeometrical cross-section to intersect D₁ and S₁. Term 3 is theattenuation along (D₁, S₁). Term 4 is the probability of interacting atthe location S₁. Term 5 is the probability of Compton scatter toward S₂.Term 6 is the solid angle for scattering toward the location S₂. Term 7is the attenuation along (S₁, S₂). Term 8 is the probability ofinteracting at the location S₂. Term 9 is the probability of Comptonscatter toward D₂. Term 10 is the solid angle for scattering toward D₂.Term 11 is the attenuation along (S₁, D₂) or (S₂, D₂). Term 12 is thedetector efficiencies.

Referring to FIG. 4, in order to properly estimate the total doublescatter contributions, we must also account for the other two distinctcontribution sources mentioned above. The estimation of the contributionfrom the emission points on the ray R²² between S₂ and D₂ is provided bythe following equation (2):

$\begin{matrix}\begin{matrix}{R_{n}^{22} = {\int\limits_{V_{S}}{\int\limits_{V_{ɛ}}{\left\lbrack {\int_{D_{2}}^{S_{2}}{{h_{n}\left\lbrack {\Delta \; {t(s)}} \right\rbrack}{n_{e}(s)}\; {ds}}} \right\rbrack_{1}\left\lbrack \frac{\sigma_{D_{2}}{dV}_{S_{2}}^{2/3}}{2\pi \; R_{S_{2}D_{2}}^{2}} \right\rbrack}_{2}}}} \\{{\left\lbrack e^{- {\int_{D_{2}}^{S_{2}}{{\mu {(E_{0,s})}}{ds}}}} \right\rbrack_{3}\left\lbrack {{\mu \left( {E_{0},S_{2}} \right)}{dV}_{S_{2}}^{1/3}} \right\rbrack}_{4}} \\{{\left\lbrack {\frac{1}{\sigma_{c}\left( E_{0} \right)}\frac{d\; {\sigma_{c}\left( {E_{0},E_{2}} \right)}}{d\; \Omega}} \right\rbrack_{5}\left\lbrack \frac{{dV}_{S_{1}}^{2/3}}{R_{S_{1}S_{2}}^{2}} \right\rbrack}_{6}} \\{{\left\lbrack e^{- {\int_{S_{2}}^{S_{1}}{{\mu {({E_{2},s})}}{ds}}}} \right\rbrack_{7}\left\lbrack {{\mu \left( {E_{2},S_{1}} \right)}{dV}_{S_{1}}^{1/3}} \right\rbrack}_{8}} \\{{\left\lbrack {\frac{1}{\sigma_{c}\left( E_{2} \right)}\frac{d\; {\sigma_{c}\left( {E_{2},E_{12}} \right)}}{d\; \Omega}} \right\rbrack_{9}\left\lbrack \frac{\sigma_{D_{1}}}{R_{S_{1}D_{1}}^{2}} \right\rbrack}_{10}} \\{{\left\lbrack e^{- {\int_{S_{1}}^{D_{1}}{{\mu {({E_{12},s})}}{ds}}}} \right\rbrack_{11}\left\lbrack {{\epsilon_{D_{2}}\left( E_{0} \right)}{\epsilon_{D_{1}}\left( E_{12} \right)}} \right\rbrack}_{12}}\end{matrix} & (2)\end{matrix}$

The estimation of the contribution from the emission points on the rayR¹² between S₁ and S₂, where the photons will each be scattered once, isprovided by the following equation (3):

$\begin{matrix}\begin{matrix}{R_{n}^{12} = {\int\limits_{V_{S}}{\int\limits_{V_{ɛ}}{\left\lbrack {\int_{D_{1}}^{S_{2}}{{h_{n}\left\lbrack {\Delta \; {t(s)}} \right\rbrack}{n_{e}(s)}\; {ds}}} \right\rbrack_{1}\left\lbrack \frac{{dV}_{S_{1}}^{2/3}{dV}_{S_{2}}^{2/3}}{2\pi \; R_{S_{1}S_{2}}^{2}} \right\rbrack}_{2}}}} \\{{\left\lbrack e^{- {\int_{S_{1}}^{S_{2}}{{\mu {(E_{0,s})}}{ds}}}} \right\rbrack_{3}\left\lbrack {{\mu \left( {E_{0},S_{1}} \right)}{dV}_{S_{1}}^{1/3}} \right\rbrack}_{4}} \\{{\left\lbrack {\frac{1}{\sigma_{c}\left( E_{0} \right)}\frac{d\; {\sigma_{c}\left( {E_{0},E_{1}} \right)}}{d\; \Omega}} \right\rbrack_{5}\left\lbrack \frac{\sigma_{D_{1}}}{R_{S_{1}D_{1}}^{2}} \right\rbrack}_{6}} \\{{\left\lbrack e^{- {\int_{S_{1}}^{D_{1}}{{\mu {({E_{1},s})}}{ds}}}} \right\rbrack_{7}\left\lbrack {{\mu \left( {E_{0},S_{2}} \right)}{dV}_{S_{2}}^{1/3}} \right\rbrack}_{8}} \\{{\left\lbrack {\frac{1}{\sigma_{c}\left( E_{0} \right)}\frac{d\; {\sigma_{c}\left( {E_{0},E_{2}} \right)}}{d\; \Omega}} \right\rbrack_{9}\left\lbrack \frac{\sigma_{D_{2}}}{R_{S_{2}D_{2}}^{2}} \right\rbrack}_{10}} \\{{\left\lbrack e^{- {\int_{S_{2}}^{D_{2}}{{\mu {({E_{2},s})}}{ds}}}} \right\rbrack_{11}\left\lbrack {{\epsilon_{D_{1}}\left( E_{1} \right)}{\epsilon_{D_{2}}\left( E_{2} \right)}} \right\rbrack}_{12}}\end{matrix} & (3)\end{matrix}$

The total double scatter contribution to the LOR (D₁, D₂) in a TOFoffset bin n is obtained by summing together the partial estimatesabove: R_(n)=R_(n) ¹¹+R_(n) ²²+R_(n) ¹², which has units of (number ofcounts)/cm²/sec. R_(n) is the ideal TOF response. To account for theactual time resolution of the detectors, a discrete time resolutionfunction, p_(mn), is applied. p_(mn) is the probability that an eventactually occurring in time bin n would be measured by the detectors intime bin m. The measured TOF response in the LOR is therefore modeled asR_(m)=Σ_(n) p_(mn)R_(n). The ideal and measured time bins may havedifferent structures. For example, in the ideal case there may be manysmall bins, while the measured data may have fewer, wider time bins.

The integrals in equations (1)-(3) are computed by sampling scatterpoint positions, but the 1/R² _(SS) factors in these equations suggestthe contributions could become arbitrarily large when the two scatterpoints get close together, making it difficult to sample themadequately. This is dealt with by using an approximate analyticalintegration over a small volume ΔV_(S) around S₁ with a specified radiusR_(S) when R_(SS)<R_(S). The possibility that R_(SS)=0 is allowed, inwhich case the scatter angle is chosen randomly. When the emissionpoints lie between S₁ and S₂, there are three terms in equation (3) thatinvolve the distance between S₁ and S₂, as shown in the first line inthe set of equations (4) below. One needs to estimate the average valueof these terms over all possible scatter point pairs lying in the spherearound S₁ with radius R_(S). Instead of estimating this by discretesampling, it is computed analytically using certain approximations. Thesecond line follows from the first under the approximation that theemitter density n_(e) and the linear attenuation coefficient μ do notvary significantly between S₁ and S₂. The third line is equivalent tothe second, and the fourth line is equivalent to the third. The fifthline follows from the fourth under the assumption that the photonattenuation over a distance R_(S) is small, and the desired averagecontribution is found to be 3n_(e)/R_(s). This value is used in place ofsampled values of these terms whenever R_(SS)<R_(S). In the TOF case,the TOF bins for the contributions are based on the positions of thesampled emission points as usual.

${{{\left\lbrack {\int_{S_{1}}^{S_{2}}{{n_{e}(s)}{ds}}} \right\rbrack \left\lbrack \frac{1}{R_{S_{1}S_{2}}^{2}} \right\rbrack}\left\lbrack e^{- {\int_{S_{1}}^{S_{2}}{{\mu {({E_{0},s})}}{ds}}}} \right\rbrack} \approx {\frac{1}{\Delta \; V_{s}}{\overset{R_{s}}{\int\limits_{0}}{\overset{\pi}{\int\limits_{0}}{\overset{2_{\pi}}{\int\limits_{0}}{\left( {n_{e}R} \right)\left( \frac{1}{R^{2}} \right)\left( e^{{- \mu}\; R} \right)R^{2}\mspace{11mu} \sin \mspace{11mu} \theta \mspace{11mu} {dR}\mspace{11mu} d\; \theta \mspace{11mu} d\; \varphi}}}}}} = {{\frac{4\pi \; n_{e}}{\Delta \; V_{s}}{\int_{0}^{R_{s}}{{Re}^{{- \mu}\; R}{dR}}}} = {{{\frac{4\pi \; n_{e}}{\Delta \; V_{s}}{\frac{1}{{\mu \left( E_{0} \right)}^{2}}\left\lbrack {1 - {e^{{- {\mu {(E_{0})}}}R_{s}}\left( {1 + {\mu \; \left( E_{0} \right)R_{s}}} \right)}} \right\rbrack}} \approx \frac{4\pi \; n_{e}R_{s}^{2}}{\Delta \; V_{s}}} = {{\frac{3n_{e}}{R_{s}}\mspace{14mu} {for}\mspace{14mu} \mu \; R_{s}}1}}}$

When the emission points lie between D₁ and S₁ or D₂ and S₂, there areonly two terms in equations (1) or (2) that involve the distance betweenS₁ and S₂, as shown in the first line in the set of equations (5) below.Following the same process as above, the desired average contribution isfound to be 3/R_(s) ². This value is used in place of sampled values ofthese terms whenever R_(SS)<R_(S).

$\begin{matrix}{{{\left\lbrack \frac{1}{R_{S_{1}S_{2}}^{2}} \right\rbrack \left\lbrack e^{- {\int_{S_{1}}^{S_{2}}{{\mu {({E_{1},s})}}{ds}}}} \right\rbrack} \approx {\frac{1}{\Delta \; V_{s}}{\overset{R_{s}}{\int\limits_{0}}{\overset{\pi}{\int\limits_{0}}{\overset{2_{\pi}}{\int\limits_{0}}{\left( \frac{1}{R^{2}} \right)\left( e^{{- \mu}\; R} \right)R^{2}\mspace{11mu} \sin \mspace{11mu} \theta \mspace{11mu} {dR}\mspace{11mu} d\; \theta \mspace{11mu} d\; \varphi}}}}}} = {{\frac{4\pi}{\Delta \; V_{s}}{\int_{0}^{R_{s}}{e^{{- \mu}\; R}{dR}}}} = {{{\frac{{4\pi}\;}{\Delta \; V_{s}}{\frac{1}{\mu \left( E_{1} \right)}\left\lbrack {1 - e^{{- {\mu {(E_{1})}}}R_{s}}} \right)}} \approx \frac{4\pi \; R_{s}}{\Delta \; V_{s}}} = {{\frac{3}{R_{s}^{2}}\mspace{14mu} {for}\mspace{14mu} \mu \; R_{s}}1}}}} & (5)\end{matrix}$

The integrals in equations (1)-(3) are computed by discrete samplingtechniques. In particular, the volume integrals over V_(S1) and V_(S2)are discretized so that each volume element dV_(S) is associated withone scatter point. For efficiency, the same set of scatter points usedfor single scatter estimation may be used for both double scatterpoints. In principle, for each dV_(S1), all dV_(S2) should be sampled.Therefore, for N scatter points, one should in principle computecontributions from all N² ordered pairs of scatter points (allowing forthe possibility of 2 scatters within one dV_(S) cell). In practice,however, complete sampling may not be practical or necessary. The pairsonly need to adequately sample the different materials present, andscattering angles. Referring to FIG. 5, it can be seen that for a firstscattering point S₁ some of the possible second scattering points, forexample S₂ and S′₂, will give very similar contributions, while others,for example S₂ and S″₂, will give very different contributions. It istherefore advantageous for the efficiency of the simulation to choose asubset of the possible scatter point pairs that reduces the redundancyof the sampling without introducing significant bias.

To this end, a reduced sampling algorithm is used: Suppose there are Nscatter points. Each scatter point has a randomly chosen position withina cell of a defined regular three-dimensional grid covering the imagingvolume of interest. There is one scatter point per cell. Let M (1≤M≤N)be a pair multiplicity such that each scatter point is to be associatedwith M other points to form M pairs per original point. The total numberof scatter point pairs to be defined is thus NM. Define a super-gridhaving approximately √{square root over (NM)} super-cells. Eachsuper-cell contains approximately N/√{square root over (NM)}=√{squareroot over (N/M)} scatter points per grid cell. Each point belongs to oneand only one cell. These √{square root over (NM)} cells are chosen to becontiguous and grouped so as to have minimal spatial extent. There areapproximately M×√{square root over (N/M)}=√{square root over (NM)} pairsoriginating in each of the √{square root over (NM)} super-cells. Foreach super-cell, the second points per pair are chosen randomly subjectto the constraint that all super-cells must be sampled, givingexhaustive sampling at the super-cell level. The scatter contributionscomputed using this reduced sampling are scaled up to account for theunder-sampling. When M→N, the super-cell is equal to a single scattercell, and the sampling is exhaustive over all possible N² ordered pairs.

FIG. 6 shows an example of the scatter grid (dotted and solid lines) andsuper-grid (solid lines only), for the case N=16, and M=4. There are 16scatter point cells and 8 super-cells, each containing 2 points. In FIG.15, one of the super-cells is shaded, and 8 possible scatter point pairsassociated with the two points in the super-cell are indicated. Out ofthe 256 possible scatter pairs in this example, only 64 will be sampled,but these are distributed with good spatial uniformity.

The DSS technique described herein can be efficiently implemented inconjunction with the SSS technique since they can share much of the sameinfrastructure, such as the scatter point distribution, detectorsampling, detector to scatter point distances and angles, emission andattenuation ray sums, detector efficiencies, the time-of-flightresolution model, photon cross-section tables, and other components. Tothese common components, the DSS technique adds the scatter pairdefinition, ray sums between the scatter points, cross-sections for thesecond scatter, analytical approximations for contributions from closelyspaced pairs, an efficient technique for sampling scatter pairs, andother components. Three new parameters are introduced to control thebehavior of the DSS: a minimum energy for tracking a doubly scatteredphoton (e.g., 350 keV), the pair multiplicity M (e.g., 100), and thescatter point separation below which the analytical approximation isused (e.g., 2.5 cm).

DSS (like SSS) uses an “absolute scaling” technique to make itsestimates of scatter quantitatively compatible with measured PET data.“Absolute scaling” means that a theoretically complete physical modelfor estimating scatter events/sec in a detector pair is used, asdescribed in equations (1)-(3). But using this to directly estimatemeasured data would require a very accurate model of the detectorefficiencies and count losses. However, the corresponding expression forunscattered events (below) involves similar efficiencies (last twoterms):

$T = {{{\left\lbrack {\int_{D_{1}}^{D_{2}}{{n_{e}(s)}{ds}}} \right\rbrack \left\lbrack e^{- {\int_{D_{1}}^{D_{2}}{{\mu {({E_{0},s})}}{ds}}}} \right\rbrack}\left\lbrack \frac{\sigma_{D_{1}D_{2}}^{2}}{2\pi \; R_{D_{1}D_{2}}^{2}} \right\rbrack}{\epsilon_{D_{1}D_{2}}^{2}\left( E_{0} \right)}}$

The ratio of these efficiencies can be estimated much more accuratelythan the individual efficiencies, so these true coincidence efficienciesare estimated using the same physics models used for the scattered eventefficiency calculations, and include in the SSS and DSS estimates:

${{e^{- {\int_{D_{1}}^{D_{2}}{{\mu {({E_{0},s})}}{ds}}}}{\int_{D_{1}}^{D_{2}}{{n_{c}(s)}{ds}}}} \sim {\left( {\left\lbrack \frac{\sigma_{D_{1}D_{2}}^{2}}{2\pi \; R_{D_{1}D_{2}}^{2}} \right\rbrack {\epsilon_{D_{1}D_{2}}^{2}\left( E_{0} \right)}} \right)^{- 1}\left( {S_{1} + S_{2}} \right)}} = {{SSS} + {DSS}}$

In other words, an estimated trues normalization is applied. This meansthat the sinogram computed by SSS+DSS is directly quantitatively(“absolutely”) comparable to a normalized trues sinogram that representsthe attenuated forward projection of the emission-rate density, withunits of events/cm²/sec.

An initial validation of the DSS algorithm was performed by comparing itto Monte Carlo (MC) simulation of 20 and 30 cm diameter cylindricalwater-filled phantoms, 26 cm long, uniformly activated with thepositron-emitting isotope ¹⁸F. These MC simulations were performed withPenelope-2011 using a cylindrical geometry approximating the Siemens mMRPET/MR scanner (detector ring diameter 65.6 cm, axial FOV 26 cm). Thecalculations assumed an energy resolution of 13.5% at 511 keV and aphotopeak energy window of 430-610 keV. These first results are for thenon-TOF case. The DSS (and SSS) simulations are absolutely scaled—thereis no cross-scaling to the MC results.

FIG. 8 shows profiles for the MC, SSS and DSS sinogram estimates for the20 cm phantom. These are radial profiles averaged over azimuthal andpolar angles, and over z. The data labeled ‘MC Trues’ are unscatteredevents, and the curve labeled ‘Model Trues’ is a simple attenuatedforward projection of the phantom. The data labeled ‘MC Tot Scat’includes all scatter orders. The curve labeled ‘SSS+DSS’ is the resultof the SSS and DSS simulations combined. The data labeled ‘MC Single’includes only single scatter events. The curve labeled ‘SSS’ is theresult of the SSS simulation. The data labeled ‘MC Double’ includes onlydouble scatter events. The curve labeled ‘DSS’ is the result of the DSSsimulation. FIG. 9 shows similar profiles for the 30 cm phantomsimulation.

Table A reports the scatter fraction estimates resulting from thevarious simulations. The scatter fraction is defined as the ratio of thesum of the specified scattered events to the sum of all scattered andunscattered events from the MC simulation. Here S₁, S₂ and S_(tot) arethe MC single, double and total scatter fractions, respectively, and SSSand DSS refer to the single and double scatter fractions from the SSSand DSS algorithms.

These results show that the SSS and DSS simulations agree very well withthe MC simulations, both in shape and in the relative magnitude ofdouble and single scatter. The scatter fractions listed in Table Aconfirm that a single plus double scatter simulation can account for allbut a few percent of the total scatter, suggesting that absolutelyscaled SSS+DSS simulations could provide a viable clinical alternativeto current rescaling algorithms.

TABLE A Scatter fractions 20 cm 30 cm MC S₁ 25.5% 30.9% SSS 25.1 ± .24%30.7 ± .12% MC S₂ 4.05% 6.96% DSS 4.05 ± .20% 6.80 ± .17% MC S_(tot)30.0% 38.9% MC (S₁ + S₂) 29.6% 37.9% SSS + DSS 29.1% 37.4% MC (S₁ +S₂)/S_(tot) 0.99 0.97 (SSS + DSS)/S_(tot) 0.97 0.96

A similar validation of DSS versus MC was performed for simulations thatincluded time-of-flight (TOF) discrimination. In this case, thesimulated phantom 100, as shown in cross-section in FIG. 17 consisted ofa 30 cm diameter and 26 cm long water-filled cylinder, in which a 14 cmdiameter by 26 cm long region, offset by 7.5 cm from the center, wasactivated with ¹⁸F. The same detector ring diameter, and photopeakenergy window, were used as previously, but in this case the energyresolution was set to 10%. A TOF resolution of 0.1 ns was assumed and 11TOF bins each 0.2 ns wide were modeled.

FIG. 18A compares the MC double scatter TOF sinograms to thecorresponding DSS TOF sinograms, for the 8 TOF bins with significantdata. In FIG. 18A, the row (A) are the MC double scatter sinograms androw (B) are the DSS TOF sinograms. FIG. 18B compares the MC totalscatter (all orders) sinograms to the corresponding SSS+DSS TOFsinograms, for the 8 TOF bins with significant data. In FIG. 18B, therow (A) are the MC total scatter sinograms and row (B) are the SSS+DSSsinograms. These sinograms represent averages over polar angle and zposition. Table B gives the ratio of the sums over these sinograms, aswell as the ratio of their total over all TOF bins (non-TOF). Thenumbers in parentheses are the standard deviations of the ratios over 10repeat SSS and DSS simulations. The good absolute agreement between theMC and SSS+DSS simulations confirms that the TOF SSS and DSS algorithmsgive an adequate representation of the physics of scattering in TOF PET,and further show that absolute SSS+DSS only slightly underestimates thetotal scatter.

TABLE B TOP offset (ns) DSS/MC double SSS + DSS/MC total −0.4 0.873(0.036) 0.853 (0.011) −0.2 1.008 (0.024) 0.971 (0.007) 0.0 1.014 (0.019)0.985 (0.005) 0.2 1.045 (0.018) 1.003 (0.006) 0.4 1.038 (0.017) 0.999(0.005) 0.6 1.027 (0.014) 0.982 (0.004) 0.8 1.018 (0.020) 0.953 (0.007)1.0 0.729 (0.031) 0.689 (0.014) non-TOF 1.023 (0.015) 0.983 (0.005)

Several tests of the DSS algorithm on clinical data were performed forthe non-TOF case, to evaluate its effectiveness and robustness. FIGS.10A-10D are coronal images from four different reconstructions ofclinical PET data using the radiotracer ¹⁸F-DCFPyL. They show that theabsolutely scaled SSS+DSS technique results in images that are arguablysuperior visually and quantitatively to the other reconstructionoptions.

FIGS. 11A-11C are zoomed image sections from three of thereconstructions shown in FIGS. 10A, 10C, and 10D, respectively, showingthat the SSS+DSS Abs algorithm gives intermediate results for scatterfractions compared to the two algorithms currently available in clinicalreconstruction software.

FIGS. 12A-12D are coronal images from four different reconstructions ofa second clinical PET study using the radiotracer ¹⁸F-DCFPyL, showingthat the SSS+DSS Abs images are comparable to those produced by theother reconstruction options.

FIGS. 13A-13C are coronal images from three different reconstructions ofa clinical PET study using the radiotracer ¹⁸F-FDG. It is found thatincreasing the double scatter pair multiplicity factor M from 100 to 500has very little effect on the image. This indicates the large values ofM, which entail longer reconstruction times, will not be necessary.

FIGS. 14A-14D are coronal images from two different reconstructions foreach of two clinical PET studies using the radiotracer ¹⁸F-FDG, showingthat in these cases SSS+DSS Abs gives comparable results to the currentstandard reconstruction SSS Rel.

FIGS. 15A-15D are coronal images from two different reconstructions foreach of two clinical PET studies using the radiotracer ⁶⁸Ga-PSMA. Theyshow that reconstructions using SSS+DSS Abs result in higher, andprobably more accurate, estimates of the scatter fraction than thecurrent clinical standard reconstruction for this type of study, SSSAbs, which is known to underestimate the total scatter.

FIGS. 16A-16B are coronal images from two different reconstructions of asecond clinical PET study using the radiotracer ⁶⁸Ga-PSMA. Again, theyshow that SSS+DSS Abs produces less underestimation of the total scatterthan the standard SSS Abs technique, because it includes double as wellas single scatter.

Applying the DSS technique described above, methods for scattercorrection in a TOF PET scanner will be described. It is to beunderstood that non-TOF DSS is included as a special case of TOF DSS.

[Double Single Scatter Estimation]—

In some embodiments, a method for applying scatter correction to a scandata acquired in a time-of-flight positron emission tomography (TOF PET)scanner assuming that both photons in the pairs of annihilation photonsare scattered once is disclosed. The method comprises:

(A) selecting a pair of detector positions in the TOF PET scanner'sdetector ring;

(B) numerically estimating the double scatter coincidence event rate ofannihilation photons in the selected pair of detector positions, whereinthe annihilation photons in a double scatter coincidence event areassumed to be scattered once, each photon by a different scatter pointin the TOF PET scanner's image volume, during a TOF PET acquisition scanperiod;

(C) selecting another pair of detector positions;

(D) repeating step (B) wherein the selected pair of detector positionsis the another pair of detector positions;

(E) repeating steps (C) and (D) until all detector pairs of interesthave been selected; and

(F) scatter correcting, or modeling scatter in, the acquired scan dataobtained from the pair of detector positions during the TOF PETacquisition scan period to reconstruct a TOF PET image based on the scandata.

Referring to equation (3), in some embodiments of the method, theselected pair of detector positions are D₁ and D₂, and the numericallyestimating step comprises:

a) selecting the two scatter points, S₁ and S₂, in the TOF PET scanner'simage volume;

b) sampling values of the emission rate density between the two scatterpoints S₁ and S₂, determining the time of flight offset bin, n, for eachcontribution, and storing these values in a vector R_(n) ¹²; [Term 1 in(3)]

c) calculating an exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₀ between thetwo scatter points S₁ and S₂; [Term 3 in (3)]

d) calculating the exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₁ between S₁and D₁; [Term 7 in (3)]

e) calculating the exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₂ between S₂and D_(z); [Term 11 in (3)]

f) calculating the reciprocals of the squares of the distances from S₁to S₂, from S₁ to D₁, and from S₂ to D₂; [Terms 2, 6 and 10 in (3)]

g) calculating the linear attenuation coefficients at energy E₀, at thepoints S₁ and S₂; [Terms 4 and 8 in (3)]

h) calculating the ratio of the differential Compton cross section for aphoton at energy E₀ to scatter to an energy E₁ to its total Comptoncross section; [Term 5 in (3)]

i) calculating the ratio of the differential Compton cross section for aphoton at energy E₀ to scatter to an energy E₂ to its total Comptoncross section; [Term 9 in (3)]

j) calculating the geometrical cross section of detector D₁ for thephoton traveling from S₁; [Term 6 in (3)]

k) calculating the geometrical cross section of detector D₂ for thephoton traveling from S₂; [Term 10 in (3)]

l) calculating the detection efficiency of detector D₁ for the photontraveling from S₁ with energy E₁; [Term 12 in (3)]

m) calculating the detection efficiency of detector D₂ for the photontraveling from S₂ with energy E₂; [Term 12 in (3)]

n) calculating the image sample volumes associated with S₁ and with S₂;[Terms 2, 4 and 8 in (3)]

o) multiplying the quantities computed in b) through n) together to getone contribution to the double scatter count rate for the detector pairD₁, D₂, for each sampled time of flight offset bin in R_(n) ¹²;

p) selecting another pair of distinct scatter points S₁ and S₂,repeating steps a) through o) and adding the results to the previousresults in R_(n) ¹²;

q) repeating step p) until the entire image volume has been adequatelysampled by both S₁ and S₂.

In some embodiments of the method, the computed line integrals, Comptoncross sections, geometrical cross sections, detector efficiencies, andscatter point positions can be cached for reuse.

In some embodiments of the method, each scatter point is chosen randomlywithin a cell of a regular spatial grid, one point per cell. Here,‘regular spatial grid’ refers to a three-dimensional rectilinear gridhaving uniform cell size. In some embodiments, a subset of all possiblescatter point pairs is used for the calculation.

[Single Double Scatter Estimation Case]—

In some embodiments, a method for applying scatter correction to a scandata acquired in a time-of-flight positron emission tomography (TOF PET)scanner assuming that one of the two photons in the pairs ofannihilation photons is scattered twice and the other of the two photonsis not scattered is disclosed. The method comprises:

(AA) selecting a pair of detector positions in the TOF PET scanner'sdetector ring;

(BB) numerically estimating the double scatter coincidence event rate ofannihilation photons in the selected pair of detector positions, whereinone of the pair of annihilation photons in a double scatter coincidenceevent is assumed not to be scattered and the other of the two photons isscattered twice, once at each of two scatter points in the TOF PETscanner's image volume, during a TOF PET acquisition scan period;

(CC) selecting another pair of detector positions;

(DD) repeating step (BB) wherein the pair of detector positions is theanother pair of detector positions;

(EE) repeating steps (CC) and (DD) until all detector pairs of interesthave been selected; and

(FF) scatter correcting, or modeling scatter in, the acquired scan dataobtained from the pair of detector positions during the TOF PETacquisition scan period to reconstruct a TOF PET image based on the scandata.

Referring to equation (1), in some embodiments of the method, theselected pair of detector positions are D₁ and D₂ and the step (BB)comprises:

aa) selecting the two scatter points S₁ and S₂ in the TOF PET scanner'simage volume;

bb) sampling values of the emission rate density between D₁ and S₁,determining the time of flight offset bin, n, for each sample, andstoring these values in a vector quantity R_(n) ¹¹; [Term 1 in (1) or(2)]

cc) calculating an exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₀ between D₁and S₁; [Term 3 in (1) or (2)]

dd) calculating the exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₁ between S₁and S₂; [Term 7 in (1) or (2)]

ee) calculating the exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₁₂ between S₂and D₂; [Term 11 in (1) or (2)]

ff) calculating the reciprocals of the squares of the distances from S₁to S₂, from S₁ to D₁, and from S₂ to D₂; [Terms 2, 6 and 10 in (1) or(2)]

gg) calculating the linear attenuation coefficients for energy E₀ at S₁and for E₁ at S₂; [Terms 4 and 8 in (1) or (2)]

hh) calculating the ratio of the differential Compton cross section fora photon at energy E₀ to scatter to an energy E₁ to its total Comptoncross section; [Term 5 in (1) or (2)]

ii) calculating the ratio of the differential Compton cross section fora photon at energy E₁ to scatter to an energy E₁₂ to its total Comptoncross section; [Term 9 in (1) or (2)]

jj) calculating the geometrical cross section of detector D₁ for aphoton pair produced between D₁ and S₁; [Term 6 in (1) or (2)]

kk) calculating the geometrical cross section of detector D₂ for thephoton traveling from S₂; [Term 10 in (1) or (2)]

ll) calculating the detection efficiency of detector D₁ for the photontraveling from the emission point with energy E₀; [Term 12 in (1) or(2)]

mm) calculating the detection efficiency of detector D₂ for the photontraveling from S₂ with energy E₁₂; [Term 12 in (1) or (2)]

nn) calculating the image sample volumes associated with S₁ and with S₂;[Terms 2, 4, 6 and 8 in (1) or (2)]

oo) multiplying the quantities computed in bb) through nn) together toget one contribution to the double scatter count rate for the detectorpair D₁, D₂, for each sampled time of flight offset bin in R_(n) ¹¹;

pp) selecting another pair of distinct scatter points S₁ and S₂,repeating steps aa) through oo) and add the results to the previousresults in R_(n) ¹¹; and

qq) repeating step pp) until the entire image volume has been adequatelysampled by both S₁ and S₂.

In some embodiments of the method, the computed line integrals, Comptoncross sections, geometrical cross sections, detector efficiencies, andscatter point positions can be cached for reuse. In some embodiments ofthe method, each scatter point is chosen randomly within a cell of aregular spatial grid, one point per cell. In some embodiments, a subsetof all possible scatter point pairs is used for the calculation.

Referring to equation (2) instead of (1), in some embodiments of themethod comprising the steps (AA), (BB), and (CC) described above, forthe case that the annihilation photons are emitted between D₂ and S₂,the step (BB) can comprise the steps aa) through qq) described above butwith the following substitutions of the variables: D₁ ↔D₂, S₁ ↔S₂,E₁→E₂, and R_(n) ¹¹→R_(n) ²². In some embodiments of the method, thecomputed line integrals, Compton cross sections, geographical crosssections, detector efficiencies and scatter point positions can becached for reuse. In some embodiments of the method, each scatter pointis chosen randomly within a cell of a regular spatial grid, one pointper cell. In some embodiments of the method, a subset of all possiblescatter point pairs is used for the calculation.

In some embodiments of the method, the total double scatter contributionfor each detector pair of interest is computed as: R_(n)=R_(n) ¹¹+R_(n)²²+R_(n) ¹². In some embodiments of the method, a discrete detector timeresolution function, p_(mn), is applied according toR_(m)=Σ_(n)p_(mn)R_(n), to estimate the probability that an eventactually occurring in time bin n would be measured by the detectors intime bin m.

According to another aspect of the present disclosure, a method forapplying scatter correction to a scan data acquired in a TOF PET scannercomprises (I) calculating the time of flight offset distribution of twodetected photons in double scatter coincidence events; and (II) scattercorrecting acquired scan data obtained during the TOF PET acquisitionscan period to reconstruct a TOF PET image based on the scan data. Insome embodiments, the step (I) can comprise:

a) selecting two detector positions D₁ and D₂ in the TOF PET scanner'sdetector ring, defining them so that D₁ is in the positive time offsetdirection;

b) creating a time offset array initialized to zero;

c) selecting two scatter points in the image volume, thereby defining adouble scatter path;

d) for each point sampled in the discrete line integrals of the emissionrate density along the segments of the double scatter path, computing(1) the distance d₁ along the double scatter path from the emissionpoint to D₁, (2) the distance d₂ along the double scatter path from theemission point to D₂, and (3) the contribution to the detector responsefrom the emission rate density at this emission point;

e) computing the time of flight offset as t=(d₂−d₁)/c, where c is thespeed of light;

f) adding the contribution to the double scatter count rate in (D₁, D₂)from this emission point into the bin of the time offset arraycorresponding to t;

g) repeating steps c) through f) until all scatter point pairs have beenadequately sampled;

h) recording the values of the time offset array in the double scattersinogram, in the appropriate time of flight bins for the detector pair(D₁, D₂); and

i) repeating steps a) through h) for all detector pairs represented inthe sinogram. In some embodiments of the method, the data in the timeoffset array are convolved with a time resolution function before beingrecorded.

According to another aspect of the present disclosure, a method forapplying scatter correction to a scan data acquired in a time-of-flightpositron emission tomography (TOF PET) scanner is disclosed. The methodcomprises estimating contributions from annihilation photons from doublescatter coincidence events arising from two scatter points S₁ and S₂when the two scatter points lie close together; and scatter correctingacquired scan data obtained from a pair of detector positions during theTOF PET acquisition scan period to reconstruct a TOF PET image based onthe scan data.

In some embodiments of the method, the estimating step comprises:defining a minimum distance R_(s) between the two scatter points S₁ andS₂; and when the distance between S₁ and S₂ is less than R_(s), using ananalytical calculation of the average contribution from certain terms inthe expressions for the double scatter event rates. These are terms 1, 2and 3 in equation (3); and terms 6 and 7 in equations (1) and (2).

In some embodiments of the method, for the case that both photons arescattered once, the three terms 1, 2 and 3 in equation (3) that involvethe distance between the scatter points are replaced by the quantity3n_(e)/R_(s), where n_(e) is the average emitter rate density at thescatter points.

In some embodiments of the method, for the case that one photon isscattered twice and the other is not scattered, the two terms 6 and 7 inequations (1) and (2) that involve the distance between the scatterpoints are replace by the quantity 3/R_(s) ².

According to some embodiments, a method for applying scatter correctionto a scan data acquired in a TOF PET scanner is disclosed. The methodcomprises: selecting a set of double scatter point pairs that contributeto double scattering of annihilation photons; and scatter correcting theacquired scan data to reconstruct a TOF PET image based on the scandata.

In the various embodiments of the method, the step of selecting a set ofdouble scatter point pairs that contribute to double scattering ofannihilation photons comprises:

a) defining a set of N single scatter points in the reconstructed PETimages from which the double scatter contribution is to be computed;

b) defining a pair multiplicity factor M which is between 1 and N;

c) defining a grid in the images having approximately √{square root over(NM)} cells, with approximately √{square root over (N/M)} scatter pointsper cell, each point belonging to one and only one cell; and

d) forming NM scatter point pairs by randomly choosing scatter pointssubject to the constraint that the scatter point pairs sample all cellpairs exhaustively. In some embodiments of the method, the N singlescatter points are uniformly distributed over the entire volume of theimages. In some embodiments of the method, the N single scatter pointsare randomly distributed over the entire volume of the images. In someembodiments of the method, the grid is a uniform rectilinear grid with aconstant cell size, covering the entire volume of the images.

FIG. 19 illustrates one embodiment of a nuclear imaging system 2. Thenuclear imaging system 2 includes a scanner for a PET modality 12provided in a first gantry 16 a. A patient 17 lies on a movable patientbed 18 that can be movable with respect to the first gantry 16 a. ThePET modality 12 includes a plurality of detectors 50 configured todetect an annihilation photons. FIGS. 3 and 4 show the detectors 50 intheir ring configuration.

Scan data from the PET modality 12 is stored at one or more computerdatabases 40 and processed by one or more computer processors 60 of anaccompanying computer system 30. The graphical depiction of the computersystem 30 in FIG. 19 is provided by way of illustration only, and thecomputer system 30 may include one or more separate computing devices.The scan data can be provided by the PET modality 12, the secondmodality 14, and/or may be provided as a separate data set, such as, forexample, from a memory coupled to the computer system 30. The computersystem 30 can include one or more processing electronics for processinga signal received from the detectors 50.

The methods and system described herein can be at least partiallyembodied in the form of computer-implemented processes and apparatus forpracticing those processes. The disclosed methods may also be at leastpartially embodied in the form of tangible, non-transitory machinereadable storage media encoded with computer program code. The media mayinclude, for example, RAMs, ROMs, CD-ROMs, DVD-ROMs, BD-ROMs, hard diskdrives, flash memories, or any other non-transitory machine-readablestorage medium, wherein, when the computer program code is loaded intoand executed by a computer, the computer becomes an apparatus forpracticing the method. The methods may also be at least partiallyembodied in the form of a computer into which computer program code isloaded and/or executed, such that, the computer becomes a specialpurpose computer for practicing the methods. When implemented on ageneral-purpose processor, the computer program code segments configurethe processor to create specific logic circuits. The methods mayalternatively be at least partially embodied in a digital signalprocessor formed of application specific integrated circuits forperforming the methods.

In some embodiments, at least one non-transitory computer-readablestorage medium is provided having computer-executable instructionsembodied thereon, wherein, when executed by at least one processor, thecomputer-executable instructions cause the at least one processor toperform embodiments of the methods described herein.

Although the subject matter has been described in terms of exemplaryembodiments, it is not limited thereto. Rather, the appended claimsshould be construed broadly, to include other variants and embodiments,which may be made by those skilled in the art.

What is claimed is:
 1. A method for applying scatter correction to ascan data acquired in a time-of-flight positron emission tomography (TOFPET) scanner, the method comprising: (A) selecting a pair of detectorpositions in the TOF PET scanner's detector ring; (B) numericallyestimating the double scatter coincidence event rate of annihilationphotons in the selected pair of detector positions, wherein theannihilation photons in a double scatter coincidence event are assumedto be scattered once, each photon by a different scatter point in theTOF PET scanner's image volume, during a TOF PET acquisition scanperiod; (C) selecting another pair of detector positions; (D) repeatingstep (B) wherein the selected pair of detector positions is the anotherpair of detector positions; (E) repeating steps (C) and (D) until alldetector pairs of interest have been selected; and (F) scattercorrecting the acquired scan data obtained from the pair of detectorpositions during the TOF PET acquisition scan period to reconstruct aTOF PET image based on the scan data.
 2. The method of claim 1, whereinthe selected pair of detector positions are D₁ and D₂, wherein the step(B) comprises: a) selecting the two scatter points, S₁ and S₂, in theTOF PET scanner's image volume; b) sampling values of the emission ratedensity between the two scatter points S₁ and S₂, determining the timeof flight offset bin, n, for each contribution, and storing these valuesin a vector R_(n) ¹²; c) calculating an exponential function of thenegative of the line integral of the linear attenuation coefficients atenergy E₀ between the two scatter points S₁ and S₂; d) calculating theexponential function of the negative of the line integral of the linearattenuation coefficients at energy E₁ between S₁ and D₁; e) calculatingthe exponential function of the negative of the line integral of thelinear attenuation coefficients at energy E₂ between S₂ and D₂; f)calculating the reciprocals of the squares of the distances from S₁ toS₂, from S₁ to D₁, and from S₂ to D₂; g) calculating the linearattenuation coefficients at energy E₀, at the points S₁ and S₂; h)calculating the ratio of the differential Compton cross section for aphoton at energy E₀ to scatter to an energy E₁ to its total Comptoncross section; i) calculating the ratio of the differential Comptoncross section for a photon at energy E₀ to scatter to an energy E₂ toits total Compton cross section; j) calculating the geometrical crosssection of detector D₁ for the photon traveling from S₁; k) calculatingthe geometrical cross section of detector D₂ for the photon travelingfrom S₂; l) calculating the detection efficiency of detector D₁ for thephoton traveling from S₁ with energy E₁; m) calculating the detectionefficiency of detector D₂ for the photon traveling from S₂ with energyE₂; n) calculating the image sample volumes associated with S₁ and withS₂; o) multiplying the quantities computed in b) through n) together toget one contribution to the double scatter count rate for the detectorpair D₁, D₂ for each sampled time of flight offset bin in; p) selectinganother pair of distinct scatter points S₁ and S₂, repeating steps a)through o) and adding the results to the previous results in R_(n) ¹²;and q) repeating step p) until the entire image volume has beenadequately sampled by both S₁ and S₂.
 3. The method of claim 2, whereinthe computed line integrals, Compton cross sections, geometrical crosssections, detector efficiencies and scatter point positions are cachedfor reuse.
 4. The method of claim 2, wherein each scatter point ischosen randomly within a cell of a regular spatial grid, one point percell.
 5. The method of claim 4, wherein a subset of all possible scatterpoint pairs is used for the calculation.
 6. A method for applyingscatter correction to a scan data acquired in a time-of-flight positronemission tomography (TOF PET) scanner, the method comprising: (A)selecting a pair of detector positions in the TOF PET scanner's detectorring; (B) numerically estimating the double scatter coincidence eventrate of annihilation photons in the selected pair of detector positions,wherein one of the pair of annihilation photons in a double scattercoincidence event is assumed not to be scattered and the other of thetwo photons is scattered twice by a scatter point in the TOF PETscanner's image volume, during a TOF PET acquisition scan period; (C)selecting another pair of detector positions; (D) repeating step (B)wherein the pair of detector positions is the another pair of detectorpositions; (E) repeating steps (C) and (D) until all detector pairs ofinterest have been selected; and (F) scatter correcting the acquiredscan data obtained from the pair of detector positions during the TOFPET acquisition scan period to reconstruct a TOF PET image based on thescan data.
 7. The method of claim 6, wherein the selected pair ofdetector positions are D₁ and D₂, wherein the step (B) comprises: a)selecting the two scatter points S₁ and S₂ in the TOF PET scanner'simage volume; b) sampling values of the emission rate density between D₁and S₁, determining the time of flight offset bin, n, for each sample,and storing these values in a vector quantity R_(n) ¹¹; c) calculatingan exponential function of the negative of the line integral of thelinear attenuation coefficients at energy E₀ between D₁ and S₁; d)calculating the exponential function of the negative of the lineintegral of the linear attenuation coefficients at energy E₁ between S₁and S₂; e) calculating the exponential function of the negative of theline integral of the linear attenuation coefficients at energy E₁₂between S₂ and D₂; f) calculating the reciprocals of the squares of thedistances from S₁ to S₂, from S₁ to D₁, and from S₂ to D₂; g)calculating the linear attenuation coefficients for energy E₀ at S₁ andfor E₁ at S₂; h) calculating the ratio of the differential Compton crosssection for a photon at energy E₀ to scatter to an energy E₁ to itstotal Compton cross section; i) calculating the ratio of thedifferential Compton cross section for a photon at energy E₁ to scatterto an energy E₁₂ to its total Compton cross section; j) calculating thegeometrical cross section of detector D₁ for a photon pair producedbetween D₁ and S₁; k) calculating the geometrical cross section ofdetector D₂ for the photon traveling from S₂; l) calculating thedetection efficiency of detector D₁ for the photon traveling from theemission point with energy E₀; m) calculating the detection efficiencyof detector D₂ for the photon traveling from S₂ with energy E₁₂; n)calculating the image sample volumes associated with S₁ and with S₂; o)multiplying the quantities computed in b) through n) together to get onecontribution to the double scatter count rate for the detector pair D₁,D₂ for each sampled time of flight offset bin in R_(n) ¹¹; p) selectinganother pair of distinct scatter points S₁ and S₂, repeating steps a)through o) and adding the results to the previous results in R_(n) ¹¹;and q) repeating step p) until the entire image volume has beenadequately sampled by both S₁ and S₂.
 8. The method of claim 7, whereinthe computed line integrals, Compton cross sections, geometrical crosssections, detector efficiencies, and scatter point positions are cachedfor reuse.
 9. The method of claim 7, wherein each scatter point ischosen randomly within a cell of a regular spatial grid, one point percell.
 10. The method of claim 9, wherein a subset of all possiblescatter point pairs is used for the calculation.
 11. The method of claim7, wherein the steps a) through Q) are performed with the followingsubstitutions of the variables: D₁ ↔D₂, S₁ ↔S₂, E₁→E₂, and R_(n)¹¹→R_(n) ²².
 12. The method of claim 11, wherein the computed lineintegrals, Compton cross sections, geographical cross sections, detectorefficiencies and scatter point positions are cached for reuse.
 13. Themethod of claim 11, wherein each scatter point is chosen randomly withina cell of a regular spatial grid, one point per cell.
 14. The method ofclaim 13, wherein a subset of all possible scatter point pairs is usedfor the calculation.
 15. A system for processing and reconstructingtime-of-flight positron emission tomography (TOF PET) sinogram data,comprising: a processor capable of executing instructions; and anon-transitory, machine readable storage medium encoded with programinstructions for controlling a TOF PET scanner, such that when theprocessor executes the program instructions, the processor performs amethod for applying scatter correction to a scan data acquired in theTOF PET scanner, the method comprising: (A) selecting a pair of detectorpositions in the TOF PET scanner's detector ring; (B) numericallyestimating the double scatter coincidence event rate of annihilationphotons in the selected pair of detector positions, wherein theannihilation photons in a double scatter coincidence event are assumedto be scattered once, each photon by a different scatter point in theTOF PET scanner's image volume, during a TOF PET acquisition scanperiod; and (C) selecting another pair of detector positions; (D)repeating step (B) wherein the selected pair of detector positions isthe another pair of detector positions; (E) repeating steps (C) and (D)until all detector pairs of interest have been selected; and (F) scattercorrecting the acquired scan data obtained from the pair of detectorpositions during the TOF PET acquisition scan period to reconstruct aTOF PET image based on the scan data.
 16. A non-transitory, machinereadable storage medium encoded with program instructions forcontrolling a time-of-flight positron emission tomography (TOF PET)scanner, such that when a processor executes the program instructions,the processor performs a method for applying scatter correction to ascan data acquired in the TOF PET scanner, the method comprising: (A)selecting a pair of detector positions in the TOF PET scanner's detectorring; (B) numerically estimating the double scatter coincidence eventrate of annihilation photons in the selected pair of detector positions,wherein the annihilation photons in a double scatter coincidence eventare assumed to be scattered once, each photon by a different scatterpoint in the TOF PET scanner's image volume, during a TOF PETacquisition scan period; and (C) selecting another pair of detectorpositions; (D) repeating step (B) wherein the selected pair of detectorpositions is the another pair of detector positions; (E) repeating steps(C) and (D) until all detector pairs of interest have been selected; and(F) scatter correcting the acquired scan data obtained from the pair ofdetector positions during the TOF PET acquisition scan period toreconstruct a TOF PET image based on the scan data.
 17. A system forprocessing and reconstructing time-of-flight positron emissiontomography (TOF PET) sinogram data, comprising: a processor capable ofexecuting instructions; and a non-transitory, machine readable storagemedium encoded with program instructions for controlling a TOF PETscanner, such that when the processor executes the program instructions,the processor performs a method for applying scatter correction to ascan data acquired in the TOF PET scanner, the method comprising: (A)selecting a pair of detector positions in the TOF PET scanner's detectorring; (B) numerically estimating the double scatter coincidence eventrate of annihilation photons in the selected pair of detector positions,wherein one of the pair of annihilation photons in a double scattercoincidence event is assumed not to be scattered and the other of thetwo photons is scattered twice by a scatter point in the TOF PETscanner's image volume, during a TOF PET acquisition scan period; (C)selecting another pair of detector positions; (D) repeating step (B)wherein the pair of detector positions is the another pair of detectorpositions; (E) repeating steps (C) and (D) until all detector pairs ofinterest have been selected; and (F) scatter correcting the acquiredscan data obtained from the pair of detector positions during the TOFPET acquisition scan period to reconstruct a TOF PET image based on thescan data.
 18. A non-transitory, machine readable storage medium encodedwith program instructions for controlling a time-of-flight positronemission tomography (TOF PET) scanner, such that when a processorexecutes the program instructions, the processor performs a method forapplying scatter correction to a scan data acquired in the TOF PETscanner, the method comprising: (A) selecting a pair of detectorpositions in the TOF PET scanner's detector ring; (B) numericallyestimating the double scatter coincidence event rate of annihilationphotons in the selected pair of detector positions, wherein one of thepair of annihilation photons in a double scatter coincidence event isassumed not to be scattered and the other of the two photons isscattered twice by a scatter point in the TOF PET scanner's imagevolume, during a TOF PET acquisition scan period; (C) selecting anotherpair of detector positions; (D) repeating step (B) wherein the pair ofdetector positions is the another pair of detector positions; (E)repeating steps (C) and (D) until all detector pairs of interest havebeen selected; and (F) scatter correcting the acquired scan dataobtained from the pair of detector positions during the TOF PETacquisition scan period to reconstruct a TOF PET image based on the scandata.